function cfdTimeAdvanceRK3
%--------------------------------------------------------------------------
%  Written by LiuHaHa @ NWPU, 2022 01 
%  Contact me at: liuzhikan@mail.nwpu.edu.cn
%==========================================================================
% Routine Description:
%   This function advance time using 3rd order Runge-Kutta
%--------------------------------------------------------------------------
global cfdFlow slvpara;
theNumberOfElements = cfdFlow.mesh.numberOfElements;
phyTemp = cfdFlow.Phy;
phyPre = cfdFlow.Phy;
elementVolumes = cfdFlow.mesh.elementVolumes;
cfdFlow.phyPreOutput = cfdFlow.Phy;

RKa(1,1)  = 1.0;      RKa(2,1)  = 0.0;      RKa(3,1)  = 1.0;
RKa(1,2)  = 3.0/4.0;  RKa(2,2)  = 1.0/4.0;  RKa(3,2)  = 1.0/4.0;
RKa(1,3)  = 1.0/3.0;  RKa(2,3)  = 2.0/3.0;  RKa(3,3)  = 2.0/3.0;
RKa(1,3)  = RKa(1,3)/RKa(1,2);

for subIrt=1:3
    
    cfdCalculateFluxAndEigenvalue;
    %
    cfdCalculateResidual;
    
    if (subIrt == 1)
        cfdCalculateDeltaT;
    end
    
    dt = cfdFlow.dt;
    resd = -cfdFlow.Resd;
    if (subIrt == 1)
        for idim=1:slvpara.ND+2 
            phyTemp(1:theNumberOfElements,idim) = phyPre(1:theNumberOfElements,idim)...
                - RKa(3,subIrt)*dt*resd(1:theNumberOfElements,idim)./elementVolumes(1:theNumberOfElements);
        end
    else
        for idim=1:slvpara.ND+2 
            phyPre(1:theNumberOfElements,idim) = RKa(1,subIrt)*phyPre(1:theNumberOfElements,idim);
            phyTemp(1:theNumberOfElements,idim) =  phyPre(1:theNumberOfElements,idim) + ...
                  RKa(2,subIrt)*phyRk(1:theNumberOfElements,idim) ...
                - RKa(3,subIrt)*dt*resd(1:theNumberOfElements,idim)./elementVolumes(1:theNumberOfElements);
        end
    end
    phyRk = phyTemp;
    cfdFlow.Phy = phyTemp;
    
end

cfdFlow.PhyTotalTime = cfdFlow.PhyTotalTime + cfdFlow.dt;
















